﻿#pragma once
#include "OCMain/FileUtil.h"
#include "OCMain/Log.h"
#include "json/json.h"
#include "rapidjson/document.h"
#include "OCGIS/ShpUtil.h"
#include "OCGIS/GDALCPP.h"
#include "mapbox/geojson.hpp"
#include "mapbox/geojson_impl.hpp"
#include "mapbox/geojsonvt.hpp"
#include "mapbox/util.hpp"
#include "mapbox/geojsonvt_io.hpp"
#include <mapbox/geometry.hpp>
#include <mapbox/feature.hpp>
#include <mapbox/variant.hpp>

#include "OCesium/TilingScheme.h"
#include "OCMain/HttpClient.h"
#define STB_IMAGE_IMPLEMENTATION
#include "OCOsg2Cesium/tinygltf/stb_image.h"

using namespace OCPP;
using namespace OCPP::Cesium;

using namespace rapidjson;

struct to_value {
	int operator()(mapbox::feature::null_value_t) {
		return 0;

	}

	int operator()(bool t) {
		return t;

	}

	int operator()(int64_t t) {
		return t;

	}

	int operator()(uint64_t t) {
		return t;

	}

	int operator()(double t) {
		return t;

	}

	int operator()(const std::string& t) {
		return 0;

	}

	int operator()(const std::vector<mapbox::feature::value>& array) {
		return 0;
	}

	int operator()(const std::unordered_map<std::string, mapbox::feature::value>& map) {
		return 0;
	}
};

inline float calcIceTerrain(float value, float B, float A)
{
	const static double e = exp(1.0);
	const static double E1 = 0.1494;
	const static double E2 = 0.3119;
	return value * (E1 * log(B) - E2) / (E1 * log(A) - E2);
}

using GridPoint = std::tuple<float, float, float>;
using GridPointVec = std::vector<GridPoint>;

FloatVector getTerrain(int z, int y, int x)
{
	long lHTTPSCode;
	HttpClientEx httpClient([](const String& message) {});
	HttpClient::UserAgent = "Mozilla/5.0 (Windows NT 6.3; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/48.0.2564.48 Safari/537.36";
	HttpClient::Referer = "www.tianditu.gov.cn";
	httpClient.SetTimeout(10);

	long retCode = 0;
	String url = str_format("http://10.104.207.155:9080/MapServer/terrarium/%d/%d/%d", z, y, x);
	String data = httpClient.fetchContent(url, retCode);

	FloatVector fVec;
	if (!data.empty())
	{
		int width, height, channels;
		auto image = stbi_load_from_memory((const uint8*)data.data(), data.size(), &width, &height, &channels, 0);
		if (image)
		{
			const static float interval = 1.0;
			const static float offset = -32768;

			fVec.resize(width * height);
			for (auto r = 0; r < height; r++)
			{
				for (auto c = 0; c < width; c++)
				{
					auto i = r * width + c;
					auto i0 = (height - r - 1) * width + c;
					double r1 = (double)image[i * 3];
					double g1 = (double)image[i * 3 + 1];
					double b1 = (double)image[i * 3 + 2];
					double height = (r1 * 256.0f + g1 + b1 / 256.0f) * interval + offset;
					fVec[i0] = height;
				}
			}
			stbi_image_free(image);
		}
	}

	return fVec;
}

void  MK_Algorithm()
{
	String fileName = "\\\\10.104.207.155\\Cesium\\Assets\\2024\\Hubei\\Shoudianqu\\24022220.021.json";
	String output = "e:/test/";

	String content = FileUtil::getFileText(fileName);
	Document document;
	if (document.Parse(content.c_str()).HasParseError()) {
		std::cerr << "Error parsing JSON!" << std::endl;
		return;
	}

	int numLon = document["numLon"].GetInt();
	int numLat = document["numLat"].GetInt();

	float startLon = document["startLon"].GetDouble();
	float startLat = document["startLat"].GetDouble();
	float endLon = document["endLon"].GetDouble();
	float endLat = document["endLat"].GetDouble();

	float dx = (endLon - startLon) / numLon;
	float dy = (endLat - startLat) / numLat;

	mapbox::geojson::feature_collection features;

	float minValue = Math::POS_INFINITY;
	float maxValue = Math::NEG_INFINITY;
	auto ds = document["DS"].GetArray();

	FloatVector gridValues;
	gridValues.resize(ds.Size());
	auto getGridValueByLonLat = [&startLon, &startLat, &dx, &dy, &numLon, &numLat, &gridValues](float lon, float lat)
		{
			float value = 0.0f;
			int c = (lon - startLon) / dx;
			int r = (lat - startLat) / dy;
			if (c >= 0 && c < numLon && r >= 0 && r < numLat)
			{
				int i = (numLat - r - 1) * numLon + c;
				value = gridValues[i];
			}
			return value;
		};
	for (auto i = 0; i < ds.Size(); i++)
	{
		auto& e = ds[i];
		float gridValue = e.GetFloat();
		minValue = std::min(gridValue, minValue);
		maxValue = std::max(gridValue, maxValue);

		int r = Math::Floor(i / numLon);
		int c = i % numLon;
		float lon = startLon + c * dx;
		float lat = startLat + (numLat - r - 1) * dy;

		if (gridValue > 0.0)
		{
			mapbox::geojson::feature feature;
			feature.geometry = mapbox::geojson::point(lon, lat);
			feature.properties = mapbox::feature::property_map();
			feature.properties["id"] = int(i);
			feature.properties["value"] = float(gridValue);
			feature.id = features.size();
			features.push_back(feature);
			gridValues[i] = gridValue;
		}
	}

	//FileUtil::writeFile(StringUtil::replaceAll(fileName, ".json", ".geojson")
	//	, mapbox::geojson::stringify(mapbox::geojson::geojson({ features })));

	mapbox::geojsonvt::Options options;
	options.maxZoom = 22;
	options.indexMaxZoom = 18;
	options.extent = 4096;
	options.tolerance = 20;
	options.indexMaxPoints = 10000;
	auto geojson = mapbox::geojson::geojson({ features });
	mapbox::geojsonvt::GeoJSONVT index{ geojson, options };
	std::shared_ptr<TilingScheme> tilingScheme = std::make_shared<WebMercatorTilingScheme>();
	int startX = 0, startY = 0;
	for (int z = 9; z <= 11; z++)
	{
		int numX = tilingScheme->getNumberOfXTilesAtLevel(z);
		int numY = tilingScheme->getNumberOfYTilesAtLevel(z);
		int y, x;
		startX *= 2;
		startY *= 2;
		int count = 0;
		int tmpX, tmpY;
		for (y = startY; y < numY; y++)
		{
			for (x = startX; x < numX; x++)
			{
				String path = Stringify() << output << "/" << z << "/" << y;
				String fileName = Stringify() << path << "/" << x << ".json";
				auto& tile = index.getTile(z, x, y);
				if (tile.features.size() == 0 || tile.num_points == 0) {
					continue;
				}
				count++;
				tmpX = x;
				tmpY = y;

				//if (z == 9 && y == 213 && (x == 415 || x == 415 || x == 415))
				{
					mapbox::geojson::feature_collection features2;
					//for (auto& f : tile.features)
					//{
					//	auto id = mapbox::geojson::identifier::visit(f.id, to_value{  });
					//	features2.push_back(features[id]);
					//}

					FileUtil::makeDirectory(path);

					auto terrainVec = getTerrain(z, y, x);
					if (terrainVec.empty())
					{
						OCPP_LOG_INFO << str_format("z(%d),y(%d),x(%d) terrain not found.", z, y, x);
					}
					else
					{
						int size = 64;
						int step = 256 / size;
						auto rect = tilingScheme->tileXYToRectangle(x, y, z);
						float west = Math::RadiansToDegrees(rect.west);
						float south = Math::RadiansToDegrees(rect.south);
						float dx = Math::RadiansToDegrees(Cesium::Rectangle::computeWidth(rect)) / (size - 1);
						float dy = Math::RadiansToDegrees(Cesium::Rectangle::computeHeight(rect)) / (size - 1);

						for (int r = 0; r < size - 1; r++)
						{
							for (int c = 0; c < size - 1; c++)
							{
								float lon = west + dx * c;
								float lat = south + dy * r;
								float gridValue = getGridValueByLonLat(lon, lat);
								if (gridValue > 0.0)
								{
									int it = r * step * 256 + c * step;
									float elev = terrainVec[it];

									float averElev = 0.0;
									for (int r1 = 0; r1 < step; r1++)
									{
										for (int c1 = 0; c1 < step; c1++)
										{
											int it2 = it + r1 * 256 + c1;
											averElev += terrainVec[it2];
										}
									}

									averElev /= (step * step);

									mapbox::geojson::feature feature;
									auto p1 = mapbox::geojson::point(lon, lat);
									auto p2 = mapbox::geojson::point(lon + dx, lat);
									auto p3 = mapbox::geojson::point(lon + dx, lat + dy);
									auto p4 = mapbox::geojson::point(lon, lat + dy);
									auto line = mapbox::geojson::linear_ring({ p1, p2, p3, p4 });

									feature.geometry = mapbox::geojson::polygon{ line };
									feature.properties = mapbox::feature::property_map();
									feature.properties["value"] = float(gridValue);
									float fixValue = gridValue;
									if (elev > 0.0 && averElev != 0.0f)
									{
										fixValue = float(calcIceTerrain(gridValue, elev, averElev));
										if (Math::isNaN(fixValue))
										{
											fixValue = gridValue;
										}
									}
									feature.properties["fixValue"] = float(fixValue);
									feature.properties["elev"] = float(elev);
									feature.properties["averElev"] = float(averElev);
									features2.push_back(feature);
								}
							}
						}

						//for (auto& f : features2)
						//{
						//	f.properties["elev"] = 1;
						//}
					}

					String jsonContent = mapbox::geojson::stringify(
						mapbox::geojsonvt::geojson{ features2 });
					size_t size = jsonContent.size();
					FileUtil::writeFile(fileName, jsonContent);

					OCPP_LOG_INFO << fileName << " num:" << features2.size() << "\n";
				}
			}
		}
		if (count == 1)
		{
			startX = tmpX;
			startY = tmpY;
		}
	}

#if 0 // Algorithm
	float x = calcIceTerrain(20.0, 2600.0, 1500.0);
	OCPP_LOG_INFO << str_format("x(%.2f)", x) << "\n";
#endif
	OCPP_LOG_INFO << str_format("min(%.2f),max(%.2f)", minValue, maxValue) << "\n";
}
